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Abstract 

We have developed a three-dimensional (3D) spectral hydrodynamic code to study 
vortex dynamics in rotating, shearing, stratified systems (e.g., the atmosphere of 
gas giant planets, protoplanetary disks around newly forming protostars). The time- 
independent background state is stably stratified in the vertical direction and has 
a unidirectional linear shear flow aligned with one horizontal axis. Superposed on 
this background state is an unsteady, subsonic flow that is evolved with the Euler 
equations subject to the anelastic approximation to filter acoustic phenomena. A 
Fourier-Fourier basis in a set of quasi-Lagrangian coordinates that advect with the 
background shear is used for spectral expansions in the two horizontal directions. 
For the vertical direction, two different sets of basis functions have been imple- 
mented: (1) Chebyshev polynomials on a truncated, finite domain, and (2) rational 
Chebyshev functions on an infinite domain. Use of this latter set is equivalent to 
transforming the infinite domain to a finite one with a cotangent mapping, and us- 
ing cosine and sine expansions in the mapped coordinate. The nonlinear advection 
terms are time integrated explicitly, whereas the Coriolis force, buoyancy terms, 
and pressure/enthalpy gradient are integrated semi-implicitly. We show that in- 
ternal gravity waves can be damped by adding new terms to the Euler equations. 
The code exhibits excellent parallel performance with the Message Passing Interface 
(MPI). As a demonstration of the code, we simulate the merger of two 3D vortices 
in the midplane of a protoplanetary disk. 
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1 Introduction 



Three of Jupiter's notable features are: rapid rotation (spin period of just 
under ten hours); an atmosphere striped with a large number of alternating 
zones and belts corresponding to strongly shearing east- west winds (hundreds 
of m/s); and many long-lived, coherent vortices, the most prominent being the 
Great Red Spot (GRS). Of course, these three characteristics - rotation, shear, 
and vortices - are all dynamically linked, so much so that it is often claimed 
that the presence of the first two implies the likely existence of the third 



sec 



Marcusl . ll99oL Il99.j for reviews). Protoplanetary disks (the disks of gas 



and dust in orbit around newly-forming protostars) also have rapid rotation 
and intense shear, which has inspired proposals that such disks shou l d also 
be populated with long-lived, coherent st orms (iBarge and Sommeria , 

Barranco and Marcusl. 



1995; 



2000, 



20051 ) . These vortices may play two critical roles in star and planet formation: 



(1) In cool, non- magnetized disks, vortices may transport angular momen- 
tum radially outward so that mass can continue to accrete onto the growing 
protostar, and (2) vortices are very efficient at capturing and concentrating 
dust particles, which may help in the formation of planetesimals, the basic 
"building blocks" of planets. 



Motivated by these geophysical and astrophysical problems, we have developed 
a three-dimensional (3D) spectral hydrodynamic code that employs specially 
tailored algorithms to handle the computational challenges due to rapid rota- 
tion, intense shear, and strong stratification. In subsonic flow, short- wavelength 
acoustic waves have periods that are much shorter than the characteristic 
timescale of the large-scale advective motions. In numerical simulations, the 
time-step for an explicit algorithm must be short enough to temporally re- 
solve these fast waves (i.e., the Courant-Friedrich-Lewy, or CFL, condition), 
which may be inefficient for calculating the evolution of the slower, large- 
scale flow for long integration times. One strategy is to filter sound waves 
from the fluid equations ( "sound-proofing" ) so that the time-step will be lim- 
ited by the longer advective timescale. The anelastic approximation does this 
by replacing the full continuity equation with the kinematic constraint that 
the mass flux be divergence-free. This approximation still allows for the ef- 
fects of density stratification (e.g., buoyancy in the vertical momentum equa- 
tion, pressure- volume work in the energy /entropy equation) and has been 
employed extensively in the study of deep, subsonic convection in planetary 
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atmospheres (lOgura and Phillipj.ll962l:lGoughLll969l : lBannon[ri 996 and stars 
fjQilman and GlatzmaierLll98lUMiesch et all l200Ct ) . In lBarranco et all (|20n(t . 

we re-derived the anelastic approximation in the context of protoplanetary 
disks. Stratified media support the propagation of internal gravity waves. As 
these waves travel from high density to low density regions, their amplitudes 
can grow to large values (so as to conserve energy flux). If the density con- 
trast is large, velocity and thermodynamic fluctuations can become sufficiently 
large so as to invalidate the anelastic approximation and/or violate the CFL 
condition. We have developed a technique based on "negative feedback" to 
damp these waves when they propagate into low-density gas which has very 
little inertia. 



We compute the evolution of the anelastic equations with a spectral method. 
The basic philosophy of spectral methods is to approximate any function of 
i nterest with a finite sum of basis functions multiplied by spectral coefficient s 
((Gottlieb and Orszagl . Il977l : iMarcusl . Il98fit ICamito et all Il988t iBovdl . Hxxj). 
A partial differential equation (PDE) in space and time is reduced to a cou- 
pled set of ordinary differential equations (ODE) for the time evolution of 
the spectral coefficients. The chief advantage of spectral methods over finite- 
difference methods is accuracy per degrees of freedom (e.g., number of spectral 
modes or number of grid points). In one dimension, the global error (e.g., L2 
norm) for a spectral method with N spectral modes scales as (1/N) , whereas 
for a finite-difference method with iV grid points, the error scales as (1/N) P , 
where p is the (fixed) order of the method. Thus, to get the same level of ac- 
curacy, spectral methods generally require far fewer degrees of freedom. This 
advantage is even more pronounced in 3D problems requiring high resolution. 



Because of the linear background shear, the fluid equations are not autonomous 
in the cross-stream coordinate, making it problematic to apply periodic bound- 
ary conditions in this direction. The equations can be made autonomous in 
the hori zontal directions by transforming to a set of Lagrangian shearing coor- 
dinat es ([Goldreich and Lvnden- Bell . 1196.4 iMarcus and Pressi . Il977t iRogallol 
19811 ). Features in the flow that are advected by the shear appear quasi- 
stationary in the shearing coordinates, allowing for larger time steps to be 
taken in the numerical integration. Because the background state generally 
depends on the vertical coordinate, we do not impose periodicity in the verti- 
cal direction. We have implemented the code with two different sets of vertical 
basis functions: (1) Chebyshev polynomials on a truncated, finite domain, and 
(2) rational Chebyshev functions on an infinite domain. Use of this latter set 
is equivalent to transforming the infinite domain to a finite one with a cotan- 
gent mapping, and using cosine and sine expansions in the mapped coordinate 



(Ca in et all Il984t iBov 



ig cosine i 
dl. l2000h . 



The code is pseudospectral. That is, nonlinear terms are not computed via con- 
volution, which is computationally slow; rather, they are computed on a grid 
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Fig. 1. Linear shear flow: v = —ayx. The system rotates counterclockwise around 
the z-axis (which points out of the page) at a rate Cl. The shear and vortex illustrated 
in this example are anticyclonic: a < 0. 

of collocation points, and then transformed to wavenumber space with Fast 
Fourier Transforms (FFTs). Different horizontal Fourier modes interact only 
through these nonlinear terms, and are otherwise independent. This motivated 
us to divide the domain into blocks along the two horizontal directions, but 
not along the vertical direction. That is, on a parallel machine, each processor 
works with only a subset of the horizontal Fourier wavenumbers, but with 
all of the vertical data corresponding to those horizontal wavenumbers. The 
only time in the simulation that different processors communicate with one 
another is in the computation of the FFTs needed for the nonlinear terms. We 
get excellent parallel performance with the Message Passing Interface (MPI) 
on the IBM Blue Horizon and Datastar supercomputers at the San Diego 
Supercomputer Center. 

In §2, we will review the anelastic equations with shear. In §3, we will describe 
the spectral method in detail. In §4, we will present tests of accuracy and 
performance, and preview potential applications. Finally, in §5, we will outline 
future work on these applications, as well as propose further improvements to 
these simulations. 



2 Anelastic Equations with Shear 

For clarity, we define up-front the following flow variables (vectors are repre- 
sented with bold type): velocity v, vorticity u> = V x v, pressure p, density 
p, temperature T, entropy s, and potential temperature 9. For an ideal gas: 
p = pTZT and s = C v ln(pp~ 7 ) = C p In 9, where C v is the specific heat at con- 
stant volume, C p is the specific heat at constant pressure, 1Z = C p — C v is the 
gas constant, and 7 = C p /C v is the ratio of specific heats. We work in a Carte- 
sian coordinate system (x, y, z) with corresponding unit vectors (x, y, z). The 
acceleration of gravity points in the negative z direction: a grav = —gz. The 
rotation axis is aligned with the z axis: CI = Clz. 



4 



Each flow variable is represented as the sum of a time-independent background 
component, denoted with overbars, and a time- varying component, denoted 
with tildes (e.g., v = v + v). The background velocity field is a unidirectional 
shear flow that is parallel to the x-axis and varies only with the y coordinate 
(see Fig. 1): 



v 



-ayx, 



where a is the constant shear rate. The background vorticity corresponding 
to this shear flow is u> = oz. The shear is cyclonic if a and Q have the same 
sign, and anticyclonic if they have opposite signs. The background thermody- 
namic state, which we assume depends only on the vertical coordinate z, is in 
hydrostatic balance: 



dp 

dz 



-pg- 



An important measure of the stratification is the Brunt- Vaisala frequency, or 
the frequency of buoyancy-driven oscillations in a convectively stable atmo- 
sphere: 



uj b {z) 



g ds 
C p dz 



g_d9_ 

e~Iz' 



(3) 



Ther e are many different variations of the anelastic approximation ( Ogura and Phillipsl . 
19621 : Gough , 19691 : Oilman and Glatzmaier , 1981 ; Bannon , 19961 ) . The key as- 
sumption they all share is that the flow is subsonic and the thermodynamic 
fluctuations are small relative to their background values: 



v \ 2 p p T 9 

1 rsj — rsj — r\j r^j — <<j^ 

c s J p p T 9 



(4) 



where e is the Mach number and c s is the local sound speed. In a typical 
derivation, all variables are expressed as asymptotic series expansions in pow- 
ers of the Mach number e, the expansions are substituted int o the governing 



fluid e quations, and terms are grouped by like powers of e. In iBarranco et al 
(2000), we derived a version of the anelastic approximation for a protoplane- 



tary disk with a constant te mperatu r e bac kground. Here, we present a more 
general version proposed by iBannon (|l996h . but modified to include a back- 
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ground shear: 



P 



;V ■ pV = V • V 



d\iap s 

dz 



v z , 



d 




dt 




d 


d 


dt 


dx 



e 



v = aVyX + v x (£j + 2VLz) — V/i + -=gz 



0- 



= -v-VQ - —v z . 

dz 

h = p/p + v 2 /2, 
p/p = p/p + f/T, 
6/6 = p/(pgH p )-p/p, 



(5a) 
(5b) 

(5c) 

(5d) 
(5e) 
(5f) 



where we have defined a dynamic enthalpy h, and a density scale height H p = 
| (din p/dz)~ \. The error in the anelastic approximation due to the neglected 
terms scales as 0(e 2 ). The careful reader may note that equation (5f) is not the 
"correct" linearization of t he potential tem perature; the denominator in the 
pressure term should be 7p. Bannon ( 1996| ) shows that the anelastic equations 
do not conserve energy unless one makes the replacement 7p — > pgH p . 

In the case where the background is constant temperature T = T , we can elim 



inate the potential tempera ture in favor of the gas temperature via ((Bannon 
199fit iBarranco et allEEnh : 



T 



' d__ d_\^ 
dt a ^ dx ) 



UJ 



-v-VT -^-{T + T)v 



(6a) 
(6b) 



It is useful for diagnostic purposes to identify energy-like conserved quantities 
for the anelastic system (5). We define kinetic, shear, and thermal energies: 

E K = [ pv 2 /2dV } (7a) 
Jv 

E s = / pv x v x dV, (7b) 
Jv 

e t = f c p pf (e/e) dv } (7c) 

where we integrate over the computational volume V: —L x /2 < x < L x /2, 
—Ly/2 < y < Ly/2, and — L z < z < L z . For this exercise, we assume that 
all quantities with a tilde are periodic in x and y. 3 We also assume that if 

3 We have not yet introduced the shearing coordinates, but if one did this exercise 
assuming that variables with tildes are periodic in the shearing coordinates, one 
would get the exact same result. 
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L z is finite, then v z — at z — ±L 2 , or if L z — > oo, then {; z — > sufficiently 
fast. To obtain an equation for the evolution of the kinetic energy, take the 
dot product of the momentum flux pv with the momentum equation (5b) and 
integrate over the domain V. Many of the resulting terms can be written as 
perfect divergences which will integrate to zero with periodicity in x and y 
and the vertical boundary conditions. The terms that survive the integration 
are source/sink terms for the kinetic energy. For the evolution of the shear 
energy, multiply the x-component of the momentum equation (5b) by pv x , 
and integrate over V. Here, the integration is not so simple because y times a 
tilde quantity is not necessarily periodic in y; this will lead to surface terms 
that do not cancel. For the evolution of the thermal energy, multiply the 
potential temperature equation (5c) by C p pT/6 and integrate over V. The 
resulting energy evolution equations are: 

dE K /dt = + £ 1 +£ 2 , (8a) 

dE s jdt = - £ 1 + £ 3 , (8b) 

dE T /dt= -£ 2 , (8c) 

d(E K + E s + E T )/dt= + £ 3 , (8d) 

where the energy source/sink terms are defined: 

£ i = / apv x v y dV, (9a) 
J v 

£ 2 = f p (9/9) gv z dV, (9b) 

£3= apv x v y Lydxdz. (9c) 

Jy=L y /2 

There are two source/sink terms for the kinetic energy: exchange with the 
background shear, and exchange with the thermal energy. As expected, these 
same terms appear, with opposite signs, in the evolution equations for the 
shear and thermal energies. We can define a "total" energy which is the sum 
of kinetic, shear, and thermal energies. The evolution equation for the total 
energy has only one source/sink term, which is the surface term from the shear 
energy evolution equation. This term represents the flow of energy into and 
out of the edges of the domain in the cross-stream direction. 



3 Spectral Method 



Let q(x, y, z, t) represent some arbitrary variable, which may be either a scalar 
or a component of a vector field. We write: 

N x /2 N y /2 N z 

q(x,y,z,t)n Yl E T,Urnn(t)e ik * x e ik «v<t> n (z), (10) 
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where £, m, n are integers, and {qe mn (t)} are the set of spectral coefficients. The 
Fourier- Fourier basis in x and y has wavenumbers: {k x = £Ak x , k y = mAk y }, 
where Ak x = 2n/L x , Ak y = 2ir/L y , and where L x and are the dimensions 
of the periodic box. The basis functions for the vertical direction, <p n (z), will 
be defined in Section 3.2. If q is a real function (and if the <p n {z) are real 
functions of z) then we have the reality condition: q* lmn = g_£_ mj „, where 
()* represents complex conjugation. We only need to keep track of half the 
spectral coefficients, say, for £ > 0. 

Our approach is a pseudospectral method. Products of two or more variables 
are not computed via convolution, which is computationally expensive. In- 
stead, we define a grid of collocation points: {xg = £Ax,y m = mAy,z n }, 
where Ax = L x /N x and Ay = L y /N y . The collocation points are uniformly 
spaced in x and y, but not in z (see Section 3.2). A variable can equally 
be represented by its set of spectral coefficients {qimn} ("wavenumber space" 
or "FFF-space") or by the values of the variable at the set of grid points 
{q(x£,y m , z n )} ("physical space" or "PPP-space"). We go between these two 
representations with the use of Fast Fourier Transforms (FFTs): 



F FT 

{qimn} {q(x£, y m , Z n )}. (11) 



To compute a product of two or more variables, we transform to physical space, 
multiply the variables at each collocation point, and then transform back to 
wavenumber space. The product will be contaminated with aliasing errors, 
but these will be small if the variables are well-resolved (that is, the spectral 
coefficients for high wavenumbers are sufficiently small). Finally, we mention 
that at times we will work in "mixed" space in which the horizontal direc- 
tions are in Fourier-Fourier space, but the vertical direction is untransformed: 
{qem(z n )}. We call this "FFP-space." 



3.1 Shearing coordinates 



The anelastic equations with shear (5) are autonomous in time t and the 
streamwise direction x, but are non-autonomous in the cross-stream direction 
y and the vertical direction z. It would be ideal if we could make the equations 
autonomous in y before imposing periodicity. Fortunately, this can be done 
by transforming to a set of quasi-Lagrangian coordinates that advect with the 
background shear ( Goldreich and Lvnden-Belll . 19651 : Marcus and Press) . 1977 ; 



x'=0 at t=0 



x'=0 at t=L x /]o[L 
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Fig. 2. Lagrangian shearing coordinates that advect with the shear. In this example, 
a < 0. Center box is computational domain of size (L x ,L y ). Eight periodic images 
are also shown. The dashed vertical lines correspond to lines of constant x' at time 
t = 0. At time t = t rm = L x /\a\L y , the periodic images, which are carried with 
the shearing coordinates, have realigned - this is the optimal time to re-map the 
sheared coordinates back onto a Cartesian grid. 



t' = t, 



x' = x + ayt, 



y = y 



z = z. 



(12a) 
(12b) 
(12c) 
(12d) 



The Jacobian for this coordinate transformation is unity; the volume of a 
computational cell is preserved as it is sheared. Partial derivatives in the two 
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coordinate systems are related via the chain rule: 



9 d d 

dl = dr +ay ^ (13a) 
I- 1 

^ ^ + t' ® (13 ) 
dy dy' dx' ' 

m = w (13d) 



The anelastic equations in the the shearing coordinates are: 







dx' I \ dz 



(14a) 



dv ( d \ ~ 9 

— = av y x + vx(u + 20£) - V' + V<rf-Qj j ~ h + $9*, (14b) 

- = - v .^ +yat '-y-- v „ (wo) 

where the vorticity is now: u) = (V + yat'd/dx') x v. The equations are 
autonomous in both x' and y', at the expense that the equations explicitly 
depend on time t' . 

Variables are represented by Fourier- Fourier expansions in the shearing coor- 
dinates, not in the original coordinates. From the point of view of the original 
coordinates, the computational grid is sheared with the background flow (see 
Fig. 2). If left unchecked, the computational grid will become greatly distorted 
in time: lines of constant x' will approach being parallel with lines of constant 
y'. Periodically, it will be necessary to re- map the shear ing coordinate system 
back onto a rectangular Cartesian grid (|Rogalldll98l . This should not be 



done at any arbitrary time. In general, variables will not be periodic with 
respect to the original coordinates except at those special times when the pe- 
riodic images line up again vertically. We define this to be the re-map time: 

(15) 



\a\Ly 

In physical space, the coordinates of a computational cell after re-mapping 
farm, y r m) can be found by inverting (12) and evaluating at the re-map time: 

x rm = x' - ay't rm = x' =F ^y', (16a) 

Vrm = V\ (16b) 
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z = L z cot© 




Fig. 3. Cosine mapping versus cotangent mapping. The coordinate £ is uniformly 
spaced from to n. Cosine mapping clusters collocation points around boundaries, 
whereas cotangent mapping concentrates collocation points around center of the 
domain. 

where the upper sign corresponds to o > and the lower sign corresponds to 
a < 0. Using the well-known Fourier shift theorem, the re-map process can 
also be done in FPP-space (x direction is in Fourier space, but the y and z 
directions in physical space): 

ql m (y, z, t rm ) = q e (y', z', t rm ) exp (ik x ayt rm ) 

= q e (y',z',t rm )exp(±2iri£y/L y ) , (17) 

where the top sign corresponds to o > 0, and the bottom sign corresponds to 
a < 0. 

We re-map the three components of velocity and the potential temperature, 
but not the dynamic enthalpy. In an anelastic code, the enthalpy is not a dy- 
namic variable; rather, it adjusts instantaneously to maintain the divergence- 
free condition on the momentum. The enthalpy should be recomputed from 
the Poisson-like equation after the momentum and potential temperature have 
been re-mapped (see Section 3.3). 



3. 2 Vertical basis functions 



In general, the background state may depend explicitly on the vertical coor- 
dinate z. Rather than imposing artificial periodici ty in the vertical direction , 
we em ploy two different domain mapping methods ([Canuto et al.Lll98 l lBovdl . 
2000). The first mapping that we use is a cosine mapping: 



L z cos£, < £ < Ti, 



which truncates the vertical domain to be within \z\ < L z . When we couple 
this mapping with a cosine series expansion, the resulting basis functions are 



11 



the Chebyshev polynomials: 



z \ , ^ , _i z 



T n f — J = cos(n£) = cos(ncos — ). (19) 
\L Z J L z 

That is, use of a cosine basis in the mapped coordinate £ is completely equiva- 
lent to using a Chebyshev basis in the original coordinate z. One disadvantage 
of this approach is that the collocation points (which are uniformly spaced in 
the mapped coordinate £) are clustered near the boundaries z = ±L 2 (see 
Fig. 3). Another disadvantage is the need to impose what may be unphysical 
boundary conditions at the boundaries z = ±L 2 . 

The s econd mapping we use is a cotangent mapping (|Cain et all 11984 iBovdl . 
2000h : 

z = L z cot 0<£<vr, (20) 
which allows us to treat the entire infinite domain — oo < z < oo. With this 
mapping, half of the collocation points are within \z\ < L z , and they are more 
concentrated around z = (see Fig. 3). The other half of the collocation points 
are widely spaced beyond \z\ > L z . We couple this mapping with both cosine 
and sine series: 

(z \ z 
— ) = cos(n£) = cos(ncot _1 — ), (21a) 
L z J L z 

S n (-^) = sin«) = sin^cot" 1 ^-). (21b) 



Bovdl (2000) calls these functions "rational Chebyshev functions" (although 



only the functions for even n are rational, whereas the ones for odd n are 
irrational functions). 

Why do we couple the cotangent mapping with both cosine and sine expan- 
sions, whereas with the cosine mapping, we use only cosine expansions? The 
derivative of the nth Chebyshev polynomial is (with L z = 1): 

dTJz) nsin«) ^ . . 

j+n odd i+n odd 

where «o = 1/2 and aj = 1 for j > 0. That is, the derivative of a Chebyshev 
polynomial is another Chebyshev polynomial of one lower degree. On the other 
hand, the derivatives of the rational Chebyshev functions are (with L z = 1): 



dC n (z) 
dz 

dS n (z) 
dz 



nsin(n£) sin 2 (£) = n{— sin[(n — 2)£] + 2 sin(n£) — sin[(n + 2)£]} 

n [S n - 2 (z) + 2S n {z) - S n+2 {z)} , (23a) 

— n cos(n£) sin 2 (£) = — n{— cos[(n — 2)£] + 2 cos(n£) — cos[(n+2)£]} 
-n \-C n - 2 {z) + 2C n {z) - C n+2 {z)\ . (23b) 
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The derivative of a rational Chebyshev function of the cosine kind is a rational 
Chebyshev function of the sine kind, and vice versa. In principle, one could 
represent the derivative of C n (z) with a convergent infinite series of C n (z), 
obviating the need for S n (z). However, we prefer to work with recursion rela- 
tions with finite numbers of terms, such as (23). Other recursion relations for 
the rational Chebyshev functions can be found in Appendix C. 

(2000) has an extensive discussion on the asymptotic behavior of the 
rational Chebyshev functions. We want to point out that the choice of cosine 
or sine kind in representing a given function is not determined by the parity of 
the function (i.e., whether it is an even or odd function), but on the asymptotic 
behavior of the function at infinity. The rational Chebyshev functions of the 
cosine kind approach their asymptotic values in even powers of 1/y, whereas 
those of the sine kind approach their asymptotic values in odd powers of 1/y. 

The product of two variables represented by the same kind of rational Cheby- 
shev functions results in a series of the cosine kind, whereas the product of 
two variables represented by different kinds of rational Chebyshev functions 
results in a series of the sine kind. Multiplying a variable by a function that is 
odd with respect to z causes its series expansion to switch kinds. Every term in 
a given equation should be represented by the same kind of expansion. Thus, 
when we use the rational Chebyshev functions for the vertical basis functions, 
v x , v y , uj z , 9, and h are represented with rational Chebyshev functions of the 
cosine kind, whereas v z , u x , and Cj y are represented with rational Chebyshev 
functions of the sine kind. 



3.3 Time integration algorithms 

In this section, we describe algorithms for integrating the anelastic equations 
forward in time. Let the hat symbol (e.g., q) denote a variable whose hori- 
zontal directions are in Fourier space (for clarity, we drop the subscripts £, m 
previously introduced). One may think of q as a either a column vector of the 
values of the variable at the different vertical positions (i.e., FFP-space), or a 
column vector of vertical spectral coefficients (i.e., FFF-space); which of these 
is intended will be clear from the context. 

The equations are integrated in time via the method of "fractional steps." 
The velocity is integrated in four steps: an advection step, a hyperviscosity 
step, an explicit pressure step, and an implicit pressure step. The potential 
temperature is integrated in two steps: an advection step and a hyperdiffusion 
step. We use a superscript N to denote a variable at the iV"th timestep (e.g., 
t N = NAt, v N = v(t = t N )). A fraction in the superscript denotes a variable at 
an intermediate stage of the integration, not after some fraction of a timestep 
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(e.g., v N+l/A is the velocity after the first of four intermediate steps). 

We have implemented two versions of the advection step, both based on an 
explicit second-order Adams-Bashforth method. In one method, the nonlinear 
advection terms are treated explicitly and the linear terms (Coriolis and buoy- 
ancy) are treated semi-implicitly. In the other method, the advection, Coriolis, 
and buoyancy terms are all treated explicitly, but with additional terms that 
damp internal gravity waves in low-density regions. For the pressure step, we 
use the semi-implicit second-order Crank-Nicholson method. We have con- 
firmed that there are no low order splitting errors, so that the algorithms are 
globally second-order accurate. 



3.3.1 Step 1: Advection — Semi-implicit treatment of Coriolis and buoyancy 
forces 



The Coriolis force couples the horizontal components of the velocity, resulting 
in inertial oscillations with the Coriolis frequency (modified by linear shear): 
use = \j2Vt{2Vt + a). Similarly, buoyancy and pressure-volume work couple 
the vertical velocity and the potential temperature, resulting in buoyant os- 
cillations (e.g. , internal gravity waves) with the Brunt- Vaisala frequency ujb 
flKundul .ll990: iPedloskv . 1979). Both inertial and buoyant oscillations can be 
described by the simple set of coupled ordinary differential equations: 



111 = +aiu 2 + fi, 

11 2 = —a 2 ui + / 2 , 



(24a) 
(24b) 



where a dot denotes time differentiation, a% > 0, a 2 > 0, and fx, f 2 are forcing 
constants. The exact solution for arbitrary time t is: 



ui(t) =«i(0) cos(o;i) + 1*2(0) 



a 1 



svn{ujt) 



+ ht 



sm(ujt) 



lot 



2(1 - cos(wt)) 



(uty 



u 2 {t) =w 2 (0) cos{ut) — m 1 (0) ( — j sin(o;t) 

+2 



sin (cut) 



out 



h 



a 2 t z 



2(1 - cos(cjt)) 
W 2 



(25a) 



(25b) 



where we have defined u 2 = aia 2 . 



We use the above solution as a template for a new integration scheme in which 
the forcing terms (e.g., f\, f 2 ) are not constants in time, but correspond to 
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nonlinear advection terms, which we treat with an Adams-Bashforth scheme: 



*4 



N 



v X uj) 



-{v ■ ve 



N-l 



v X uj) 



AT 



jv-r 



The velocity and potential temperature are updated via: 



Vx 4 



f 2 4 



:cos(w c At)i£ + AtqiiucAt) [(2Q + <j)?)f + 91 
+ ±At 2 g 2 (cu c At)(2fi + a) -i(k' y + atk' x )h J 
cos(uj c At)Vy + Atq^cucAt) [(-2ft)u* + ^ 
+ ±At 2 g 2 (w c At)(-2fi) " 



- ifcl/i 



cos(w B A*)uf + At qi (u; B At) (g/6)6 N +% 
+ \At 2 q 2 {uj B At){g/Q) [Wl 
§ N +\ =cos{uj B At)9 N + At qi (u B At) [{-d6/dz)v? + Wl 
+ \At 2 q 2 (uj B At)(-de/dz) [m z - dh N /dz 

where the functions q\ and q 2 are defined: 

gi(r) = sin(r)/r, 
g 2 (r) = 2[1 - cos(r)]/r 2 . 



(26a) 
(26b) 



(27a) 
(27b) 
(27c) 
(27d) 



(28a) 
(28b) 



The time differencing error for one step scales as At 3 . For uc/ B At <C 1, q\ 
and q 2 can be replaced with unity without changing how the error scales with 
the size of the step. 



3.3.2 Step 1: Advection — Explicit treatment of Coriolis and buoyancy forces 
with gravity-wave damping 

As an internal gravity wave propagates from a high density region into a 
low density region, its amplitude increases so as to conserve energy flux. If 
the density contrast is large, velocity and thermodynamic fluctuations can 
become sufficiently large so as to invalidate the anelastic approximation and/or 
violate the CFL condition. We have developed a strategy to control gravity 
waves that is akin to "negative feedback" in circuit theory. Mathematically, 
the origin of the buoyant oscillations is the fact that the time derivative of 
the vertical velocity depends linearly on the potential temperature, and vice 
versa. Physically, a hot parcel of fluid rises into less dense material, expands 
and cools; it then sinks, compresses, and heats-up, and the cycle continues. 
To damp these oscillations, we add a term proportional to the time derivative 
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of the potential temperature to the vertical velocity equation, and vice versa. 
This introduces a term that is proportional to — v z into the equation for the 
vertical velocity, and a term that is proportional to —6 into the equation for 
the potential temperature. In the following integration scheme, we add the 
Coriolis and buoyancy terms to the nonlinear advection terms and treat them 
all with an Adams-Bashforth scheme: 



2 



N 



— N Q 

[v X u>) - 2VLz x v N + av!/ x + -^-qz 
v ' y q J 



N-1 



[v X w 



2ttz x v + aijy x 



-gz 



(29a) 



971 




N-1 



(29b) 



The velocity and potential temperature are updated via: 
v N+1 ^=v N + At ™ ' W ^ 




<92 



(30a) 
(30b) 



The function (3(z) determines the range of z at which the damping operates; 
P(z) = for most of the domain, but will be nonzero and small in the low 
density regions. 

We note that one could damp buoyant oscillations just by adding a term 
proportional to — v z (instead of 86/ dt) to the vertical velocity equation, and/or 
a term proportional to —6 (instead of —8v z /8t) to the potential temperature 
equation. However, the advantage of adding time derivatives is that as the flow 
approaches a steady-state, the damping naturally turns off. We will explore 
the effects of this gravity-wave damping mechanism in Section 4.3. 



3.3.3 Step 2: Hyperviscosity 



Unlike finite-difference codes which have a "grid viscosity" associated with the 
spatial differencing, spectral codes have no inherent grid dissipation. Energy 
will cascade down (via nonlinear interactions) to the smallest resolved scales 
(highest wavenumbers) where it will accumulate and eventually cause large 
truncation and aliasing errors. To control this, artificial viscosity or some other 
kind of low-pass filtering is appli e d to dissipate the energy at the smallest 
scales f Canuto et all 1198 1 iBovdl . Il998l ). This method is related to "large- 
eddy" simulations in which the large scale motions are fully resolved, whereas 
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the small scale flow is treated with a sub-grid scale turbulence model. 



For the two horizontal Fourier directions, we apply a hyperviscosity of the 
form: v h ^ p {— 1) P+1 V^, where v h ^ v is a hyperviscosity coefficient, and p is an 
integer between 1 and 6. This is simple to implement because the hypervis- 
cosity operator is an exact eigenoperator of the Fourier basis functions: 



■ iy +1 \7 2 P exp(ik x x + iky-y) = — {k 2 + k 2 ) p exp(ik x x + ik y y). 



(31) 



The vertical direction is not resolved with a Fourier basis and is not periodic. 
Applying a real diffusion operator to the vertical direction is problematic be- 
cause the second derivative operator (or its powers) is not an eigenoperator of 
the Chebyshev polynomial or rational Chebyshev function bases. Furthermore, 
a real diffusion operator would raise the order of the differential equations and 
require additional unp hysica l vertic al boundary conditions and/or create spu- 
rious boundary layers. iBovdl (1998) advocates replacing the second derivative 
operator with an eigenoperator of the basis functions. For example, for the 
Chebyshev polynomial basis (where we assume L z = 1): 



dz 



dz 2 
d 



vT 



dz 



dz 



dz 



TJz) = -n TJz) 



(32a) 
(32b) 



Likewise, there is a similar replacement for the eigenoperator of the rational 
Chebyshev function basis. The chief advantage is that because the eigenop- 
erators are singular at the boundaries, no additional unphysical boundary 
conditions are needed and no spurious boundary layers are created. 

In practice, we apply hyperviscosity every time step: 



f hyp = expI-A^A? + u^n 2p )}, 



~N+- 
V * 

nN+1 



fhyp: 
fhypi 



(33a) 
(33b) 
(33c) 



where k\ = k' 2 + (k' y + otk' x ) 2 and n is the order of the Chebyshev polynomial 
or rational Chebyshev function. The hyperviscosity coefficients are initially 
set to values such that the highest resolved Fourier or Chebyshev number has 
an e-folding time equal to one timestep. They can be dynamically adjusted 
every few hundred timesteps so that the energy spectrum does not curl-up at 
the highest wavenumbers. 
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3.3.4 Step 3: Poisson-like equation for pressure/ enthalpy 



Subtraction of the enthalpy gradient is done with a semi-implicit Crank- 
Nicholson (globally second-order in time) algorithm: 



~ 7V+- 

V + 4 



v N+ i - \AtV N h N , 



~ N+l ~ N+ 
V = V 



t - ±AtV 



(34a) 
(34b) 



where V 



[ik' x , i(k' y + at N k' x ), d/dz]. Step (34a), the explicit step, uses 
the current enthalpy, whereas step (34b), the implicit step, uses the future 
enthalpy. We require that the final velocity satisfy the anelastic constraint, 
which yields a Poisson-like equation for the final enthalpy: 



v" +1 + 



d\ap s 

dz 



v 



N+l 



dz 1 



+ 



d\np\ d 

dz / dz 



h N+l 



= 0, 
2 

~ At 



where (A;^ +1 ) 2 = k' x 2 + (k' y + at N+l k' x f 



~N+1 

V + 



<ilnp N 

dz 



(35a) 



~ N+- 



(35b) 



If the vertical domain is finite (i.e., Chebyshev polynomial basis), then explicit 
boundary conditions will be needed to solve this Poisson-like equation. We 
require that the vertical velocity vanish at the walls. Imposing this on the 
^-component of (34b) yields: 



\N+1 



dh N+1 



dz 



z=±L z 



Z = ±L; 



= o, 



2 „jv+§ 
At 



z=±L z 



(36a) 
(36b) 



If the vertical domain is infinite (i.e., rational Chebyshev function basis), then 
no explicit boundary conditions are needed to invert the Poisson-like equation 
because the vertical basis functions individually satisfy the boundary condi- 
tions (i.e., Galerkin method with natural boundary conditions): dC n {z)/dz — > 
and S n (z) — * as z — > ±oo. 

The Poisson-like equation is singular for the case of k' x = k' y = 0. The anelastic 
condition for this mode implies that there can be no net vertical momentum 
across any horizontal plane: 



v N+1 
dz 



0, 

aT 2 ' 00 • 



(37a) 
(37b) 
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The integration constant for h N+1 is determined by the condition that the 
total mass in the domain is conserved. 



One can write the Poisson-like equation (35b) in matrix form: M.h = g, where 
M is a matrix containing the Laplacian-like differential operator, h is a column 
vector of vertical spectral coefficients for the enthalpy, and g is a column vector 
of vertical spectral coefficients for the right-hand side of (35b). When needed, 
the last two rows of the matrix M are overwritten with the boundary condition 
data (36b), at the expense that the divergence condition will not be satisfied 
for the two highest vertical spectral modes. This is usually an exponentially 
small error, but if it is undesirable, one can use the tau-method with Greens 
functions to correct t he di vergence at the two highest spectral modes (see 
Appendix B; iMarcusl Il984h . 



The speed and efficiency in inverting the Poisson-like equation depends on the 
functional form of the stratification term dhxp/dz. If the stratification term 
is particularly simple (e.g., zero for no stratification, a constant for constant 
stratification, or proportional to height z as in a protoplanetary disk), then one 
can use Chebyshev or rational Chebyshev recursion relations; the matrices will 
be banded and consequently quick to invert. In particular, for the finite domain 
with the Chebyshev polynomial basis, the matrix equation can be rewritten 
as: Ah = Bq, where A and B are pentadiagon al matrices (see Appendix A, 



Gottlieb and Orszagi . ll977tlCanuto et all 119881 ). If the stratification term is an 



odd function of z, then the even and odd Chebyshev modes decouple, reducing 
the pentadiagonal matrix equation to two tridiagonal matrix equations for 
the even and odd Chebyshev modes. For the infinite domain with rational 
Chebyshev function basis, the matrix M is nonadiagonal (see Appendix C). If 
the stratification term is an odd function of z, then the even and odd rational 
Chebyshev modes decouple, reducing the nonadiagonal matrix equation to 
two pentadiagonal matrix equations for the even and odd rational Chebyshev 
modes. On the other hand, if the stratification term is a more complicated 
function of height, then the terms multiplied by the stratification term will 
have to be computed in physical space on a grid of collocation points, and then 
transformed back to function space. The matrices in the Poisson-like equation 
will consequently be full and slower to invert. 



4 Numerical Tests 



4-1 Parallelization & timing analysis 



We have developed our simulation for use on parallel supercomputers such as 
the Blue Horizon (decommissioned in April 2004) and DataStar supercomput- 
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Fig. 4. Average wall-clock or real time to compute one timestep using 16, 32, 64, 
and 128 processors on the IBM Blue Horizon at the San Diego Supercomputer 
Center. Parallelization is implemented with the Message Passing Interface (MPI) 
protocol. 'Tot' = Total time for one timestep, 'Adv (comp)' = computation time for 
advection step, 'Adv (MPI)' = MPI communication time for advection step, 'Pres' = 
computation time for pressure step, 'Hyper' = computation time for hyperviscosity 
step. The dashed lines are power-law fits; the numbers to the right of the data are 
the best-fit exponent: Time oc N p / oc . For reference, we note that the DataStar runs 
approximately two times faster than the Blue Horizon. 



ers at the San Diego Supercomputer Center. The DataStar has 176 nodes of 
8 processors (CPU speed of 1.5 GHz and peak performance of 6.0 GFlops). 
Each node shares 16GB of memory. 

Different horizontal Fourier modes interact only through the nonlinear advec- 
tion terms, whereas the vertical spectral modes (Chebyshev polynomials or 
rational Chebyshev functions) are fully coupled at every stage of the compu- 
tation. This motivated us to divide the computational domain only along the 
horizontal directions. Each processor works with a small subset of horizontal 
data, but all of the corresponding vertical data: in physical space, the j'th pro- 
cessor would contain data for: x{ < x < x 3 2 , y{ < y < yi, —L z < z < +L Z ; in 
wavenumber space, a given processor would contain data for: k 3 x l < k x < k 3 x 2 , 
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fcj^i < ky < ky 2 , < n < N z . The only time different processors communicate 
with each other is when a horizontal FFT needs to be computed. There are 
generally two common strategies for this: one could use a FFT that works 
across processors, or one could transpose the data so that a given direction 
is not split across processors. For example, if one wanted to compute a FFT 
along the x direction, one would swap the data so that each processor had 
a subset of the y and z data, but all of the corresponding x data. We have 
implemented this second strategy with an efficient transpose algorithm devel- 
oped by Alan Wray at NASA- Ames/Stanford Center for Turbulence Research 
(Wray, personal communication). 

Figure 4 shows the results of timing tests for a version of the code that used 
the Chebyshev polynomial basis. The tests were done on the Blue Horizon 
and analyzed with the Vampir parallel trace tool. For reference, we note that 
the DataStar runs approximately two times faster than the Blue Horizon. We 
present the average wall-clock or real time to compute one timestep for 64 3 , 
128 3 , and 256 3 domains on 16, 32, 64, and 128 processors. We fit power laws 
(T oc Np r P c ) to the timing data and determined how the times for different 
steps in the algorithm scale with number of processors (ideal performance is 
p « 1). As expected, the advection step is the most time-intensive step in 
the algorithm 80-90% of the time is spent in this stage alone due to the large 
number of FFTs required to compute the nonlinear terms. The computation 
time ('Adv comp') scales well with the number of processors, whereas the 
communication time ('Adv MPF) has a much flatter scaling. 



4-2 2D vortex dynamics in a shearing flow 



The simplest way to validate the shearing box a lgorithm is to confirm that 
the code correctly simulates 2D vortex dynamics. iMoore and Saffmanl ()197lh 
analytically determined the steady state solutions for 2D elliptical vortices of 
uniform vorticity u z embedded in a uniform shear of strength a, and showed 
that the aspect ratio of a vortex, \ = A^/A^, where and A y are the major 
and minor axes of the ellipse, was determined by the ratio of the vorticity to 
the shear: 



a 



(3* 



\x- V x 

which is valid for vortices that rotate in the same sense as the shear. As 
one would expect, stronger vortice s are more round and compact, whereas 
weaker vortices are more elongated. Kida ()198ll ) showed that if such a vortex 
was perturbed, its major axis would oscillate about the x-axis, and its aspect 
ratio would oscillate about its mean value (area remaining constant). These 
motions would be damped via the occasional shedding of thin filaments of 
vorticity which would eventually be dissipated by viscosity. 
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Figure 5 shows the kinetic, shear, and total energies as a function of time 
for the evolution of a 2D elliptical vortex. The parameters for these runs 
were: (L x , L y ) = (4,2), (N x ,N y ) = (256,128), SI = 1, a = -1.5, t rm = 4/3, 
At = t rm /100, (A x ,A y ) = (1,0.25), u g (0) = -0.625. The initial vortex was 
surrounded with a weak halo of oppositely-oriented vorticity so that there 
was no net circulation around the vortex (r = j A u z dxdy = 0). The Kida 
oscillations are readily apparent in the energy plots as the larger amplitude, 
longer period oscillations (~ 8 re-map periods). As expected, increasing the 
hyperviscosity resulted in the damping of these oscillations; the final state is 
a steady-state vortex whose major axis is aligned with the x-axis and whose 
aspect ratio and area are constant. 

Superposed on the Kida oscillations is another low amplitude, short period 
oscillation, which is most clearly seen in the plots of total energy. This period 
is exactly the re-map period and is due to interactions of the vortex with its 
periodic images, which are realigned with the primary vortex every re-map 
period (see Fig. 2). The amplitude of these oscillations decreases with domain 
size. 



In Figure 6, we plot the percent error in kinetic energy after one re-map period 
as a function of timestep size for the explicit and semi-implicit treatments of 
the Coriolis (modified by shear) terms. We fit power-laws and confirmed that 
both schemes are second-order accurate in time. The semi-implicit treatment, 
however, has an error that is a order of magnitude smaller than the explicit 
scheme. 



4-3 Linear eigenmodes in a stratified protoplanetary disk 



To test the vertical dimension of the code, we investigate the behavior of 
linear eigenmodes (i.e., internal gravity waves) in a stratified background. 
The background shear makes it difficult to compute eigenmodes because the 
linearized equations depend explicitly on ayd/dx in the original coordinates, 
or explicitly on atik x in the shearing coordinates. If we restrict the following 
analysis to the case of d/dx = ik x = 0, then the process of finding eigenmodes 
is straight-forward. We look for eigenmodes of the form q(z) exp(ik y y — iut), 
where q represents any variable. The real part of the eigenvalue is related 
to the phase speed of the wave, c = TZ£(cu) /k y , whereas the imaginary part 
7 = TM.{uj) corresponds to wave growth (if positive) or damping (if negative). 



22 



14 

12- 
10 

8 

6 



x 10 



hyp 



■■ 10" 



Kinetic 
Shear 
Total 



oj 4 - 

C 
LU 

2\ 
0-' 
-2- 
-4- 



x 10 



hyp 



■■ 10" 



Kinetic 
Shear 
Total 



UJ 



aj 4 - 

2-, 
0-\ 
■2- ( , 
4- 







20 40 60 80 

Time (number of remap periods) 



100 







20 40 60 80 

Time (number of remap periods) 



Fig. 5. Kinetic, shear, and total energy for 2D elliptical vortex embedded in a 
linear shear. The larger amplitude, longer period (~ 8 re-map periods) oscillations 
are associated with the major axis of the vortex oscillating about t he a;- axis and with 
the aspect ratio oscillating around its mean value, as described by Kidal (1981)- The 
lower-amplitude, shorter period oscillations (exactly equal to the re-map period) 
are due to interactions with the periodic image vortices, which are realigned with 
the primary vortex every re-map period. The Kida oscillations are damped via the 
shedding of thin filaments which are dissipated by hyperviscosity. 



100 



The linearized equations are: 



—iuv x = (2Q + a)v. 



-iuVy = —2Qv x — ikyh, 
-iojv z = —dh/dz + (0/9)g — (3ujbVz, 
-iu6 = -(d6/dz)v z - (3uj b [6 - {6/g)dh/dz 
= ikydy + dv z / dz + (d In pj dz)v z , 



(39a) 
(39b) 
(39c) 
(39d) 
(39e) 



where we have included the damping terms described in section 3.3.2. If the 
damping terms are neglected, then this system of equations can be simplified 
to the second-order eigenvalue problem for the vertical velocity fluctuation: 



-)' 



£ ky 



C ^y 



(40) 



where we have denned the second-order linear differential operator: 



dz 



df_ + / dlnp \ j 
dz \ dz J 



(41) 
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Fig. 6. Fractional error in final kinetic energy after one re-map period versus 
timestep size for explicit and semi-implicit treatment of the Coriolis force. Both 
methods are second-order accurate, but the error is an order of magnitude smaller 
for the semi-implicit method. 

If the eigenvalue problem is solved on the finite domain with the Chebyshev 
polynomial basis, then boundary conditions must be imposed: v z (z = ±L z ) = 
0. On the infinite domain with the rational Chebyshev function basis, bound- 
ary conditions are naturally satisfied by the basis functions (i.e., Galerkin 
method). 

In order to proceed, we need to specify the background state. Here, we choose 
conditions relevant to a protoplanetary d i sk of gas and dust in orbit around a 
newly-formed protostar ( Barranco et all l2000l : iBarranco and Marcusl . 2005). 



The velocity profile of the gas is very nearly Keplerian: Vx(r) = rf^(r) 



GM+/r, where G is the gravitational constant, M* is the mass of the proto- 
star, and r is the cylindrical radius to the protostar. The horizontal shear rate 
of this base flow is: ok = r(dVLx / dr) = —(3/2)Qk- Keplerian shear is anti- 
cyclonic (as indicated by the negative sign), and is comparable in magnitude 
to the rotation rate itself. We do not simulate the entire disk; rather, we we 
simulate the hydrodynamics only within a small patch (Ar <C r , A<p <C 2n) 
that co-rotates with the gas at some fiducial radius r with angular velocity 
Q = Q(r ). We map this patch of the disk onto a Cartesian grid: r — r — > y, 
r Q (4> — 4>q) — > —x, z — > z, v r — > v y , — > —v x , and v z — > v z . In this rotating 
patch, the velocity profile of the gas is very nearly a linear shear flow with 
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(T = <J K (ro) = — (3/2)f2 . The Coriolis frequency (modified by the shear) is 

The vertical component of the protostellar gravity in this patch of the disk 
is very nearly g z (z) m Q^z. We assume that the background is constant tem- 
perature T = T , which yields a hydrostatic density profile that is Gaussian: 

p(z) = po exp(-z 2 /2H 2 ), Hi = nT /Q 2 , (42a) 
H p {z) = | (dlnp/dz)- 1 1 = Hl/\z\. (42b) 

The Brunt- Vaisala frequency is: 

co B (z) = ^n/C P n \z\/H . (43) 

The disk can be divided into two regimes: a nearly unstratified midplane 
(uj b < Q , \z\/H < 1.6) and a strongly stratified atmosphere above and 
below the midplane (uo B > \ z \/Hq > 1.6). 

Figure 7 shows the vertical velocity component of the undamped {(3 = 0) 
eigenmodes for k y Ho = n. The eigenvalues are purely real. There are two kinds 
of eigenmodes: those confined to the weakly stratified midplane of the disk with 
frequencies u < Q Q , and those confined to the strongly stratified regions near 
the boundaries with frequencies uo > f2 . The latter class of eigenmodes do 
not exist when the domain is infinite (rather, there is a continuous spectrum 
of waves whose velocity amplitudes are formally unbounded as \z\ — > oo). 

To investigate the effect of damping on these eigenmodes, we need to specify 
a profile for j3(z), for example: 

P(z) = I3 max {l + tanh[a(,2 - z c )/H ] + tanh[a(,2 c - z)/H ]}, (44) 

which is nearly zero for most of the domain, but rises sharply at \z\ m z c 
to its maximum value /3 maa; ; the parameter a sets how steep the rise is. In 
the calculations presented here, we take a = 4 and z c = 3H . Table 1 and 
Fig. 8 show the effect of varying f3 max on the two kinds of eigenmodes. Those 
eigenmodes that are localized around the midplane are only weakly damped, 
whereas the eigenmodes that are localized in the strongly stratified regions 
near the boundary are strongly damped. 

In fully nonlinear 3D calculations, this damping is necessary for numerical 
stability. As internal gravity waves propagate from a high density regions into 
low density regions, their amplitudes increase so as to conserve energy flux, 
causing velocity and thermodynamic fluctuations that can become sufficiently 
large so as to invalidate the anelastic approximation and/or violate the CFL 
condition. 
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Fig. 7. Vertical velocity component of the undamped (/3 = 0) eigenmodes in a strati- 
fied protoplanetary disk for (k x , k y ) = (0, tt/Hq). There are two kinds of eigenmodes: 
those localized around the weakly stratified midplane with frequencies to < Qq, an d 
those localized in the strongly stratified regions near the boundaries with frequencies 
to > Qq. 
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Table 1 



Effect of damping on linear eigenmodes. An eigenmode that is localized near the 
midplane (e.g., u = 0.8597) is weakly damped, whereas an eigenmode that is local- 
ized far from the midplane (e.g., to = 1.229) is strongly damped. Damping rate is 
linearly proportional to (5 m ax ■ 
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Fig. 8. Effect of damping on linear eigenmodes. Eigenmodes that are localized 
near the midplane (lower group) are weakly damped, whereas eigenmodes that are 
localized far from the midplane (upper group) are strongly damped. Damping rate 
is linearly proportional to f} max . 

4-4 Merger of two 3D vortices in a protoplanetary disk 



As a preview of potential applications, we present the simulation of a merger 
of two 3D vortices in the midplane of a stratified protoplanetary disk. A more 
complet e discussion of s i mulat i ons of vortices in protoplanet ary disks can be 
found in B arranco et al. (l2nnnh:lBarranco and Marcusl (IgOOfih . We use the 2D 
vortex solution of Moore and Saffmanl (|l97l| ) and lKida (|l981 ) to initialize an 
elliptical vortex in the midplane at z — 0. The vorticity is extended off the 
midplane according to a Gaussian profile: 



uj z (z 



(45) 



We initialize the temperature so that the buoyancy force exactly balances the 
vertical pressure force: 

f = I^ d MA (46 ) 
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Of course, only under very special circumstances would this procedure just so 
happen to also exactly balance the temperature equation. In general, dT/dt ^ 
initially, leading to an immediate evolution of the temperature. Vertical mo- 
tions will then be generated by the temperature changes through the buoyancy 
force, and these vertical velocities will then couple the horizontal motions at 
different heights. We allow such a 3D vortex to evolve and relax to a quasi- 
equilibrium. We then use that vortex as a template to construct an initial 
condition with two vortices destined to merge. 

Fig. 9 shows an example of a merger. The first column shows x — y slices 
in the midplane z = of the z-component of vorticity. Blue corresponds to 
anticyclonic vorticity, red corresponds to cyclonic vorticity. The second column 
shows isovorticity surfaces for the z component of vorticity, in blue, and vortex 
lines (lines that are everywhere tangent to the vorticity vector) in red. The 
vorticity field is divergence-free (since it is the curl of another vector field). 
Thus, like magnetic field lines, vortex lines cannot begin or end in free space, 
they must either end on a boundary, extend to infinity, form closed loops, or 
form unending open curves. 



5 Conclusion 

We have developed a three-dimensional, spectral, hydrodynamic code to study 
vortex dynamics in rapidly rotating, intensely sheared, and strongly stratified 
systems such as giant planet atmospheres and protoplanetary disks. The code 
integrates a number of specially tailored algorithms to handle the numerical 
challenges associated with shear and stratification. We used the anelastic ap- 
proximation to filter sound waves and shocks. The horizontal coordinates are 
transformed to a set of shearing coordinates that advect with the background 
shear flow, and a Fourier-Fourier basis in these shearing coordinates is used for 
spectral expansions in the two horizontal directions. For the vertical direction, 
two different sets of basis functions have been implemented: (1) Chebyshev 
polynomials on a truncated, finite domain, and (2) rational Chebyshev func- 
tions on an infinite domain. Use of this latter set is equivalent to transforming 
the infinite domain to a finite one with a cotangent mapping, and using co- 
sine and sine expansions on the mapped coordinate. The nonlinear advection 
terms are time integrated explicitly, whereas the enthalpy gradient and the 
terms responsible for internal gravity waves are integrated semi-implicitly. We 
show that gravity waves can also be damped by adding new terms to the 
Euler equations. The code is parallelized with the Message Passing Interface 
(MPI), and we get excellent parallel performance on the IBM Blue Horizon 
and Datastar supercomputers at the San Diego Supercomputer Center. 

The Chebyshev polynomial code has the disadvantage that it clustered the 
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Fig. 9. The merger of two 3D vortices in the midplane of a protoplanetary disk. 
In the first column, we show x — y slices in the midplane z = of the z-component 
of vorticity. Blue corresponds to anticyclonic vorticity, red corresponds to cyclonic 
vorticity. In the second column, we present isovorticity surfaces for the z-component 
of vorticity in blue, and vortex lines (lines that are everywhere tangent to the 
vorticity vector) in red. Units of time: r or & = 2tt/Qq = 1 yr. 
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collocation points near unphysical walls. This motivated us to implement the 
the rational Chebyshev functions for the infinite domain. This removed the 
walls from the problem, but it did not significantly improve the spacing of 
the collocation points. Although the collocation points were clustered more 
towards the center, only half of all the collocation points were in the "active 
region" \z\ < L z , whereas the other half were wasted on resolving the regions 
toward infinity \z\ > L z , where our equations are no longer valid. Also, the 
rational Chebyshev function expansions converged poorly if any fluid dynamics 
occurred near the mapping parameter L z , forcing us to choose a larger value 
of L z . We could improve the performance of the rational Chebyshev code if we 
used momentum ariables instead of velocity variables, since the multiplication 
by the mean density made the momentum decay faster with height. 
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to thank Andrew Szeri, Xylar Asay-Davis, and Sushil Shetty useful comments 
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A Solving the Poisson-like equation for the enthalpy with Cheby- 
shev polynomial basis ("tau method") 



For the case where the stratification is a linear function of height z, the 
Poisson-like equation for the enthalpy (35) can be written in the form: 



Mf(z) 



— - z— - A 

dz 2 ^ dz 

dj_ 

dz 



z=±l 



(A.la) 
(A.lb) 



where we have chosen dimensions so that the boundaries are located at z — 
±1. There will be one such equation for each pair of horizontal wavenumbers 
{k' x ,k' y }. We expand f(z) and g(z) in Chebyshev polynomial series: 



N N 

f( z ) = Yl fnTn(z) and g(z) = 9nT n (z). 



(A.2) 



?1=0 



n=0 
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One can then use Chebyshev recursion relations to write the above Poisson-like 
equation in matrix form: 



Mf = g, 



N 



(A.3a) 
(A.3b) 



n=l 



where M is a (A+l) x (JV+1) matrix for the second-order linear differential 
operator A4, and / and g are column vectors of N+l Chebyshev spectral coeffi- 
cients. There are two problems with this formulation: the matrix M is singular 
and the system is overspecified: we have N+l undetermined spectral coeffi- 
cients, but N + 3 constraints (N+l equations +2 boundary conditions). The 
tau method fixes both problems simultaneously by discarding the equations 
for the two highest Chebyshev modes and replacing them with the boundary 
conditions, resulting in the differential equation not being satisfied for those 
same two highest modes. The matrix M is upper triangular (except for the 
last two rows which now contain boundary conditions) and fast to invert, but 
unfortunately is not diagonally dominant. The largest elements occur down 
the last c olumn, and inverting it in its c urrent form wi l l have large round- 
off errors. Gottlieb and Orszae (1977) and Canuto et al. ( 1988j ) show how to 
rewrite the Poisson-like matrix equation in the form: 



Af = Eg, 



(A.4) 



where A an d B are diagonally do minant, pentadiagonal matrices. Gottlieb and Orszad 



( 19771 ) and lCanuto et al. 



give the recursion relations for the case 7 = 0. 



Here, we present the recursion relation for the more general case: 

C n -2fn-2 + Cnfn + Cn+2fn+2 = d>n-29ri-2 + d n g n + c/ n+ 2^ n +2, 

for 2 < n < A, 



(A.5) 



where we have defined: 

a n _ 2 A + (n- 2)7 



C n -2 '■ 

c n = 1 + 

C n +2 '- 

and 



4n(n — 1) 

/? ra+2 (A- 7 )-(l-/3n+2)7(rc + l) 
2(n 2 - 1) 

(3 n+i \ - (3 n+2 {2 - p n+4 h(n + 2) 

An(n + 1) 



2 for n = 0, 
1 for n 7^ 0, 



and (3 n 



dn-2 
dn+2 



OLn-2 



An(n — 1) ' 

fin+2 

~2(n 2 - 1)' 
4n(n + 1)' 



1 for n < N, 
for n > N. 



(A.6) 



(A.7) 
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The boundary conditions can be rewritten: 

N _|_ N _ 

£ »7n = and £ n*f n = *+-*^. (A.8) 

n=2 A n=l 2 

n even n odd 

Note that the even and odd Chebyshev polynomials are decoupled in the 
matrix equation and in the boundary conditions, so one could in fact write 
the pentadiagonal system as two sets of tridiagonal systems for the even and 
odd modes. 



B Greens functions method to correct divergence condition for two 
highest Chebyshev modes 



In Appendix A, we presented a method to solve the Poisson-like equation for 
the enthalpy in a Chebyshev polynomial basis, subject to the constraint that 
the vertical velocity vanish at the two boundaries at z = ±1. The last two rows 
of the Poisson-like matrix equation were overwritten with boundary condition 
data, resulting in the anelastic divergence condition not being satisfied for 
the two highest Chebyshev modes. Here, we present a method to co rrect the 
diver gence condition at the two highest modes using Greens functions ( Marcusl . 
1984), at the expense that equation for the evolution of the vertical momentum 
will not be satisfied at those same two highest modes. 

To make the presentation more clear, we will first define a number of operators: 



D = (d/dz), (B.la) 

D A = (d/dz) + (dlnp/dz), (B.lb) 

V* +1 = ik' x x + i(k' y + at N+1 k' x )y + Dz, (B.lc) 

7 A N+1 = ik' x x + i{k' y + <rt N+l k' x )y + D A z, (B.ld) 

Ar = v/ +1 .V M . (B.le) 



In the fourth fractional step (34b), we subtract the implicit half of the enthalpy 
gradient as before, but we also add two additional functions, each consisting 
of exactly one Chebyshev mode: 

. N+1 = q N+ i _ ^+1^+1 + ^iv+i^ + r iV+i rM ) ) (B .2) 

where we have defined n ee (At/2)h, and where T\ and T2 are scalars (as usual, 
we suppress the subscripts for wavenumber dependence). To find the enthalpy 
at the A^ + l step, we impose the anelastic constraint Va ■ v N+1 = 0, which 
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yields: 



JV+l 



V 



* + t? +1 D a T m ^ + r? +i D A T M , 



N+l 



DYi N+ \z = ±l) = vr*(z = ±l) + +1 Tm_i(±1) + rf +1 T M (±l). 
We can break this up into three separate Poisson-like equations: 



n 



7V+1 



frN+l 



"T '1 1 1 T '2 1 2 



JV+lpiV+1 



where 1^ and T 2 are Greens functions, and where: 



^A ii 



(z = ±l) 

A N+l r N+l 
A 1 1 

DT? +1 {z = ±l) 

A N+1 T N + 1 

DT% +1 (z = ±l) 



v A N+1 .v N +i : 
^ +f (,=±i), 

DaTm-i, 
Tm-i(±1), 
DaTm, 
T M (±1). 



(B.3a) 
(B.3b) 

(B.4) 



(B.5a) 

(B.5b) 
(B.5c) 
(B.5d) 
(B.5e) 
(B.5f) 



Each of these 3 Poisson-like equations can be solved via the tau method as 
in Appendix A. By construction, Va • v N+1 equals zero for all Chebyshev 
modes except the two highest; for those modes, the values of the anelastic di- 
vergence are functions of rf +l and r^ v+1 . We use these two degrees of freedom 
to impose that the anelastic divergence equal zero at the two highest modes 
as well, which yields the 2x2 matrix equation: 




where we have defined: 
a n 

«12 
«21 
"22 

A = [V A 



V T 2 



A N + 1 T N+1 _ 
A N + 1 T N + 1 _ DaTm 
A N + l T N + l _ DaTm _ x 

A? +1 r? +1 - D A T M 



(B.6) 



A 1 2 



n=M-V 
n=M-l' 
n=M' 
n=M' 



V 



^A ii 



n=M-l' 



n=M 7 



(B.7a) 
(B.7b) 

(B.7c) 
(B.7d) 

(B.7e) 
(B.7f) 



where the notation q represents the Mth Chebyshev mode of function 



n=M 



q. Note that if M is even, then T^ +1 is an even polynomial, is an odd 

polynomial; if M is odd, then r^ +1 is odd and r^ +1 is even. Using these 
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symmetries, one can show that an = a 22 = 0, and t x +1 = (3 2 /a 2 i, r 2 +1 = 
Pi/a 12 . 

This method effectively triples the computational time for the enthalpy step 
(three Poisson-like solves instead of one per timestep). However, the enthalpy 
step is still only of order 10% of the total computational time, as the FFTs 
in the nonlinear advection dominate the computation. We also note that if we 
were not working in Lagrangian coordinates (which introduced explicit time 
dependence into the gradient operators), then the Greens functions r\ and Y 2 
(two for each set of Fourier wavenumbers) would be independent of time, and 
would only have to be computed once in a pre-processing step, and stored. 



C Recursion relations for rational Chebyshev functions 



Here we present some useful recursion relations for the rational Chebyshev 
functions. We define: 



N N 



m^EfnCn(z) Or f(z) = J2fnS n (z), 
n=0 n=0 
N N 

g(z) = J2^C n (z) or g(z) = J2 g s n S n (z). (C.l) 



n=0 



n=0 



The first derivative operator converts a C n (z) series into a S n (z) series, and 
vice versa. The recursion relations have the same form for either case, except 
for the lowest two modes. Let g = df/dz, then: 



qS/C 



-(n-2)fZH + 2nfW s - (n + 2)/„% 5 ] /(4L Z ), for n > 2, 
^ = 0, g? =(3j?-3/f) /(4L Z ), OR 

= f 2 s /(2L z ), g? = (-/f + 3f 3 s ) /(4L Z ). (C.2) 



The second derivative operator preserves the kind of expansion; that is, the 
second derivative of a C n (z) series is another C n (z) series, and the second 
derivative of a S n (z) series is another S n (z) series. The recursion relations for 
the second derivative operator have the same form for either case, except for 
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the lowest four modes. Let g = d 2 f /dz 2 , then: 



[n 



+ 4(n 



-4)(n 
2)(n4 



2)/^? + 4(n 



So* 
9? 

9 S s 



8f 2 C 



8/f 



-6/p + 21# - 15£ 



9/? 
0, 



2)(n-l)^f 
n + 4)(n + 2)/2J 

/(16^), 
/(16^), 

54/ 3 c + 80/f - 35/ 7 c l /(16L 2 Z ), OR 



/(16^). 
c 



?c 



-24/ 2 c + 48/f-24/ 6 



c 



-6/f + 27/f 



15/ 5 5 



-24/ 2 5 + 48/f-24/ 6 s 
7/f - 54/f + 80/ 5 5 - 



/mi), 

35/fl /(16L*). 



for n > 4, 



(C.3) 



Like the second derivative operator, the z(d/dz) operator preserves the kind 
of expansion. Let g = z(df/dz), then: 



fJ S = 
ft 



( n -2)fZ! S 2 +(n + 2)til° 2 /4 for n > 2, 



■c/s 



/f/2, 
0, 



#1 



c 



/f + 3/ 3 c ] /4, 
■A s + 3//I /4, 



OR 



(C.4) 



One last important operation is multiplication by z. A simple recursion rela- 
tion exists only for multiplication by z on a S n (z) series, which turns it into a 
C s (z) series. Let g(z) = zf(z), then: 



( N 



£ /jf for n = 0, 

k even 

JV 

/„ 5 + 2 E fk f o r « even, 

k even 
TV 

/ n 5 + 2 ]T fk S fornodd. 



(C.5) 



k=n 
k odd 



In practice, the recursion relations for z{d/dz) and d 2 /dz 2 are unnecessary, as 
matrices for these operations can be built up from the first derivative matrices 
and the multiplication by z matrix. 
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